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Abstract. We discuss the kinematical structure of a two-parameter family of isotropic models with a central black hole. The 
family contains the slope of the central density cusp and the relative black hole mass as parameters. Most of the basic kinematical 
quantities of these models can be expressed analytically. This family contains three distinct models where also the distribution 
function, differential energy distribution and spatial LOSVDs can be expressed completely analytically. Each of these models 
show a drastically different behaviour of the distribution function. Although the effect of a black hole on the distribution function 
is very strong, in particular for models with a shallow density cusp, the differential energy distribution is only marginally 
affected. We discuss in detail the effects of a central black hole on the LOSVDs. The projected velocity dispersion increases 
with black hole mass at small projected radii, but the effect of a black hole on the shape of the LOSVDs (characterized by 
the h 4 parameter) is less straightforward to interpret. Too much reliance on the wings of the LOSVDs and the value of the h 4 
parameter to determine black hole masses might hence be dangerous. 
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1. Introduction 

During the past decade, various observational discoveries have 
changed our view on galactic nuclei rather drastically. Firstly, 
galactic nuclei are generally observed to have a cuspy den- 
sity distribution at small radii, with densities behaving as 
p(r) oc r~ r . The observed cusp slopes are far from uni- 
form and range from zero to more than 2 (Lauer et al. 1995; 
Gebhardt et al. 1996; Ravindranath et al. 2001; Seigar et 
al. 2002; Genzel et al. 2003; Scarlata et al. 2004). Secondly, 
high-resolution imaging has revealed that a substantial frac- 
tion of the galactic nuclei in both spiral and elliptical galaxies 
show small-scale structure, in the form of bars, mini-spirals or 
dust lanes (Phillips et al. 1996; Malkan, Gorjian & Tam 1998; 
Carollo, Stiavelli & Mack 1998; Tran et al. 2001; Martini 
et al. 2003). Thirdly, there is now enough credible evidence 
that (nearly) all nearby galaxies harbor a supermassive black 
hole in their centre. Intriguingly, the masses of these puta- 
tive black holes are tightly coupled to large-scale parameters 
of the host galaxies (Kormendy & Richstone 1995; Ferrarese 
& Merritt 2000; Gebhardt et al. 2000; Graham et al. 2001; 
McLure & Dunlop 2002; Ferrarese 2002; Baes et al. 2003a; 
Marconi & Hunt 2003). Clearly, these various aspects are not 
isolated features, they all play part in the processes which shape 
galactic nuclei and galaxies in general. Unfortunately, we still 
know little about the formation processes and evolutionary sce- 
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narios of galactic nuclei and their black holes. For example, it is 
still unclear when and how supermassive black hole are being 
fuelled and what the mutual influence is between black holes 
and the central density cusps in galactic nuclei. Various relevant 
processes probably play at scales which are beyond (or at the 
limit of) the current observational capabilities. In order to in- 
crease our understanding, more detailed theoretical modelling, 
for example using N-body of SPH simulations, is necessary. 

As a starting point or general framework for such studies, 
one needs a set of reference models which are simple enough 
and still can present a wide enough variety in structural char- 
acteristics. Because of the observed cuspy nature of galaxy 
centers, well explored models such as the Plummer model 
(Plummer 1911; Dejonghe 1987) and the isochrone sphere 
(Henon 1959, 1960) are less suitable. Scale free models are 
a useful alternative. Such models have a density profile that de- 
creases as r y , and they include the singular isothermal sphere 
as a special case. Due to the simplicity of these models, their 
dynamical structure can be easily studied, even in axisymmet- 
ric or triaxial generalizations (e.g. Toomre 1982; Evans 1994; 
Qian et al. 1995; de Bruijne et al. 1996; Evans et al. 1997; Jalali 
& de Zeeuw 2001). Unfortunately, scale free models always 
have an infinite total mass, with the mass diverging in the cen- 
tre for y > 3 and at large radii for y < 3. 

A set of models that does not suffer from this disad- 
vantage is the family of y-models, introduced independently 
by Dehnen (1993) and Tremaine et al. (1994). These mod- 
els have a central r~ y density cusp (y < 3) and the den- 
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sity falls at r~ 4 at large radii, such that the total mass is al- 
ways finite. Special cases are the well-known Hernquist and 
Jaffe models (Hernquist 1990; Jaffe 1984). Many interesting 
dynamical properties such as the intrinsic and projected ve- 
locity dispersions, the distribution function and the differen- 
tial energy distribution can be calculated analytically for the y- 
models, under the assumption of a self-consistent isotropic dy- 
namical structure. Extensions towards an anisotropic distribu- 
tion functions or flattened geometry have been presented (e.g. 
Dehnen & Gerhard 1994; Hiotelis 1994; Carollo et al. 1995; 
Zhenglu 2000; Baes et al. 2002a; Zhenglu & Moss 2002). 

As supermassive black holes appear to the present in nearly 
all galactic nuclei, it would be very interesting to extend the 
family of y-models with a central black hole. In principle this is 
quite straightforward: one just has to add an extra contribution 
from the black hole to the stellar potential and re-calculate the 
dynamical properties with this new potential. In fact, Tremaine 
et al. (1994) consider the case where a central black hole is 
present in the y-models. They present analytical expressions 
for the velocity dispersions and discuss the effect of a black 
hole on the distribution function. Unfortunately, they do not 
seek for analytical expressions for the distribution function and 
they leave some of the most interesting kinematical properties 
such as the differential energy distribution and the LOSVDs 
undiscussed. Ciotti (1996) provides a first attempt to construct 
completely analytical dynamical models for galaxies with a 
central black hole. Extending the work by Carollo et al. (1995), 
he considers a set of Hernquist models embedded in a dark mat- 
ter halo. By setting the dark halo radius to zero this halo reduces 
to a central black hole. He demonstrates that many of the in- 
teresting dynamical quantities, including the distribution func- 
tion and the differential energy distribution, can be calculated 
analytically (albeit as rather complicated functions involving 
Jacobian functions). Baes & Dejonghe (2004) present a one- 
parameter family with a steep y = | cusp slope and the black 
hole mass as a parameter. Almost all interesting properties of 
this family of models can be written in terms of elementary 
functions. 

In this paper, we present a detailed analysis of a two- 
parameter family of spherical isotropic dynamical models 
based on the y-models. The family contains as parameters the 
slope of the central density cusp y and the ratio // of the central 
black hole mass to the total mass of the system. The param- 
eter space covered by these models goes from weakly cusped 
models to very centrally concentrated models with an infinite 
stellar potential well, and from self-consistent models without 
black hole to systems where the dynamical structure is com- 
pletely dominated by the central black hole. In Section 2 we 
define the models. In Section 3 we derive some basic proper- 
ties, most of which can be calculated completely analytically. 
In Section 4 we discuss the energy budget of the models and 
look at the virial theorem. In Section 5 we derive expressions 
for and discuss the distribution function and the differential en- 
ergy distribution. In Section 6 we study the LOSVDs of the 
models in our family and discuss the observational signature of 
a black hole. In the last Section we summarize the results, and 
in the Appendices we present some mathematical expressions 



which might be useful for people who wish to use these models 
as input for further theoretical studies. 

2. Definition of the models 

The y-models have a luminosity density 
3-y 1 



47T 
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The parameter y determines the density slope of the system at 
small radii; it can assume values between and 3. All models 
have a luminosity density that behaves as r~ 4 at large radii, such 
that the total luminosity is finite. 

The gravitational potential of the models we consider is the 
sum of two contributions: the stellar mass and a central black 
hole. We introduce the parameter ,u as the relative importance 
of the black hole mass to the total mass in the system, such 
that ji can assume values only between and 1 . Note that the 
convention we use is similar to the one adopted by Baes & 
Dejonghe (2004), but is different from the convention used in 
e.g. Tremaine et al. (1994) and Zhao (1996). In these papers, 
li denotes the black hole mass relative to the stellar mass and 
the normalization is such that the stellar mass is set to unity. We 
prefer to set the total mass of the galaxy equal to unity however, 
because all models then have the same behaviour at large radii. 

In the limit yi — 0, there is no black hole and we recover the 
self-consistent y-models described in detail by Dehnen (1993) 
and Tremaine et al. (1994). The potential reduces to the stellar 
potential 
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This expression is not valid for y = 2, the Jaffe model. If we 
take the limit y — > 2 in the previous expression we obtain the 
Jaffe potential *P t (r) = ln(l + 1/r). This special case divides 
the family of y-models in two classes: self-consistent y-models 
with y < 2 have a finite potential well with depth = 1/(2 - 
y), whereas the stellar potential well is infinitely deep for the 
models with 2 < y < 3. 

The other extreme case on the range of possible black hole 
masses is fi = 1, corresponding to systems where the entire 
mass resides within the central black hole. In this case, the total 
gravitational potential reduces to a Kepler potential, 
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Although these systems where the dynamics are dominated by 
a central black hole do not form a realistic representation of real 
galaxies, they are useful in order to investigate the maximal 
effect of a black hole on the dynamical properties of the y- 
models. 

The general case < /i < 1 is intermediate between these 
two extreme cases. We can write the cumulative mass function 
and the potential of the y-models as 

M(r) = (l-/z)( T ^-) 3 7 +n, (4) 
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The two equations Q and Q form a two-parameter family 
of potential-density pairs. Each pair of parameters (y,p) com- 
pletely defines a dynamical model. The study of these models 
as a function of the parameters y and p is the goal of this paper. 

3. Basic properties 

3.1. The surface brightness 

The surface brightness profile of a y-model can be found by 
projecting the luminosity density on the plane of the sky, 



/(*) = 2 



f 

Jr 



p(r) r dr 
Vr 2 - R 2 ' 



At large radii, the surface brightness falls as R 3 . In the cen- 
tral regions, models with y < 1 have a finite central surface 
brightness (in spite of the divergence of the luminosity den- 
sity), whereas models with y > 1 have a surface brightness pro- 
file that diverges as R l ~ y . The surface brightness can be eval- 
uated analytically for all integer and half-integer values of y. 
Analytical expressions and other photometric quantities such 
as the cumulative surface brightness and the half-light radius 
can be found in Hernquist (1990), Dehnen (1993), Tremaine et 
al. (1994) and Baes & Dejonghe (2004). 

3.2. The velocity dispersion 

For a spherical isotropic system, the intrinsic velocity disper- 
sion profile can be found using the solution of the lowest-order 
Jeans equation (Dejonghe 1986; Binney & Tremaine 1987), 



<T 2 (r) 



1 r 

~p(r)J r 



M(r) p(r) dr 
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After substitution of the general cumulative mass function @ 
and some algebra, we obtain the expression 

p(r)cr 2 (r) = (1 - fiW y {r) + 2p W m {r), (8) 

' 2 

where the function W y (r) is defined as 

w = irl oT^- (9) 

From equation © we see that W 7 (r) is nothing but the velocity 
dispersion p(r)cr 2 (r) of the self-consistent y-model, which can 
be evaluated in terms of elementary functions for all values of 
y (see Appendix A). 

3.3. The projected velocity dispersion 

The projected velocity dispersion profile cr p (R) can be found 
by projecting the intrinsic dispersion on the plane of the sky, 



2 _ j_ C°° p(r)cr 2 (r)rdr 
P ~ W) Jr 
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This expression can be written in a form very similar to equa- 
tion M, 



I(R)ot(R) = (1 - fi)Y r (R) + 2ju Y^iR), 
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with (Tremaine et al. 1994) 
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The function Y y (R) represents the projected velocity dispersion 
I{R)cr 2 J/) of the self-consistent y-model. It can be expressed 
analytically for all integer and half-integer values of y, and has 
to be calculated numerically for the other values of y. A useful 
transformation for the numerical calculation of the integral can 
be found in Appendix B of Dehnen (1993). 



( 6 ) 4. The energy budget and the virial theorem 



The (scalar) virial theorem states that any steady-state system 
satisfies the relation 

2K=U, (13) 

where K is the total kinetic energy and U = -W is the total 
(binding) potential energy of the system. It is easy to verify 
that the y-models without a central black hole satisfy the virial 
theorem. For a self-consistent spherical system the total poten- 
tial and kinetic energy can be found as 



U, =2n p(r) ¥,(>-) r 2 dr. 



p(r)crl(r) r 2 dr 



K, = 6n f 
Jo 



2tt p(r)M,(r)rdr, 
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(14) 
(15) 
(16) 



where the second expression for the kinetic energy is derived 
by substituting equation into the first expression and par- 
tial integration. Through substitution of the expressions ([0, (0 
and (|4} into these formulae, one finds that the self-consistent 
y-models have an infinite energy budget for y > |, whereas for 
y < 4 we obtain (Tremaine et al. 1994) 



U, = 2K* 



1 
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For systems with a central black hole, the total kinetic energy 
can still be calculated using equation J15i or (I16> . with the con- 
tribution of the black hole taken into account in the dispersion 
or the cumulative mass distribution. The total potential energy, 
however, cannot be found by just replacing the stellar potential 
by the total potential, as the expression (I14t is derived under 
the assumption that the system is self-consistent. The correct 
expression for the potential energy consists of an internal and 
an external contribution (Binney & Tremaine 1987, problem 
8.2), 
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^DT P (r)Xi, »« (r)dr + fff P (r) k-VYextMldr. (18) 



In our present case, the external potential is the gravitational 
potential from the central black hole, such that these formulae 
reduce to 



U = 2tt(1 



-A0 p(r)*F, 
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(r)r dr + Anp p(r)rdr, (19) 
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or equivalently, 



where t = f( v F») is defined as 



Jo 



U = 2n\ p(r) r 2 dr + 2np. \ p(r)rdr. 
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Comparing this expression to the equivalent formula HAl of the 
self-consistent models, we notice that the black hole mass adds 
an extra contribution to the potential energy apart from its con- 
tribution to the potential of the system. Using the expressions 
0, Q and one finds that the total energy budget is infinite 
for y > 2, whereas for y < 2 we obtain 



U = 2K ■■ 



1 



1 



5-2y 

The virial theorem is thus satisfied. 
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5. The distribution function and the differential 
energy distribution 

All the kinematical information on a given system is contained 
in the distribution function /(r, v), which represents the number 
density of stars in six-dimensional phase space. For spherical 
isotropic systems, the distribution function depends only on the 
binding energy £ = *P(r) - v 2 /2. The key to calculating the 
distribution function /(£) of isotropic spherical models is the 
augmented density p(T'), i.e. the luminosity density written as a 
function of the potential. The Eddington formula specifies how 
the distribution function can be calculated from the augmented 
density, 

d_ r e dp 

dfiJo d¥ 
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d^ 2 
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The distribution function obviously is an important character- 
istic of the dynamical structure of galaxies. However, it is not 
straightforward to physically interpret the meaning of the dis- 
tribution function when expressed as a function of binding en- 
ergy. A more natural diagnostic quantity is the differential en- 
ergy distribution (DED) N{&), which describes the number of 
stars per unit binding energy. For isotropic systems, the DED 
is simply the product of the distribution function /(£) and the 
density of states function g{&), defined as the phase space vol- 
ume accessible for a star with binding energy fi. This function 
can be calculated as 



/-To 

g(8) = 16 V2tt 2 J 
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dT 
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and depends only on the potential of the system, not on the 
density profile (Binney & Tremaine 1987). 

5.1. Models without a central black hole (ju = 0) 

For y-models without a black hole, we can calculate the distri- 
bution function /*(£) directly by substituting the equations 
and into the Eddington relation i22\ . We obtain the expres- 



sion 
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This integral can generally be expressed in terms of hypergeo- 
metric functions, 

MS) = V5[(4 - y) iFi (l, ^; |; (2 - y)S) 

- 2(3 - y)2fi(l,^;|; (2 -y)S) 

+ 2(l-y) 2 F 1 (l,g;|;(2-y)£) 

+ y2 F 1 (l,g;|;(2-y)£)]. (27) 

For y — 2 ± 1 fn with n a natural number, the distribution func- 
tion can be written as a combination of elementary functions. 
In particular, the self-consistent y-models with y — 1, | and 
| have fairly simple distribution functions (Hernquist 1990; 
Dehnen 1993; Tremaine et al. 1994; Baes & Dejonghe 2004). 
The y — 2 model is a particular case: its distribution function 
can be expressed most conveniently in terms of the error func- 
tion and Dawson's integral (Jaffe 1983). 

The calculation of g(S) is also straightforward for the self- 
consistent y-models. Substituting the potential in equa- 
tion (1241 one immediately obtains (Dehnen 1993) 



(^0 

g.(S) = \b<2i? y 



(28) 



with t = f( v F«) as in equation i26i . The integral in this equation 
is quite similar to the integral in equation i25\ . A closed expres- 
sion for g*(£) for general y in terms of hypergeometric func- 
tions cannot be obtained, but the integral can be expressed ana- 
lytically for all y — 2 + 1 jn with n a natural number. Examples 
of such closed expressions can be found in Dehnen (1993) and 
Baes & Dejonghe (2004). 

5.2. Models dominated by a central black hole (jj. = 1) 

For models where the potential is completely dominated by the 
central black hole, the distribution function can also be calcu- 
lated by a straightforward application of Eddington's formula. 
The result can generally be expressed in terms of hypergeomet- 
ric functions, 

MS) = 4 1 g " y) £ 5/2 [63 2 F, (3, 6 - y; 2; -&) 

-36(1 -y)G 2 F l (3,6-y;|;-£) 

- 4y(l - y) fiVi (3, 6 -y; -£;-£)]. (29) 

For half-integer or integer values of y, the distribution function 
can be written in terms of elementary functions. In the central 
regions of the galaxy (£ — > oo), this distribution function has 
the asymptotic behaviour 



f(8) 3 " y F(r+1) ^-3/2 
Mb > 2(2n)^T(y-l) ' 



(30) 
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It is interesting to compare these results with those obtained by 
de Bruijne et al. (1996). They consider a set of axisymmetric 
cuspy densities in a spherical potential and present two differ- 
ent analytical families of three-integral distribution functions. 
For the special case of isotropic spherical galaxies dominated 
by a black hole potential (q - 6 - I and p 1 — 0) the expres- 
sion J30i is recovered, after correction for the different normal- 
ization conventions. 

Since the density of states function only depends on the 
potential of the system, we recover the simple and well-known 
expression g,{&) — V2~7r 3 £~ 5 ^ 2 , independent of y. 

A particularly interesting case is the model with y - \- 
This model has a simple distribution function and differential 
energy distribution, 



/.(£) 
N.(8) = 



V8 S 5 ' 2 



7T 3 (1 +«5) 5 ' 
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(31) 
(32) 



For -/-models with a central black hole, this model is the one 
with the smallest possible cusp slope. Indeed, y-models with 
y < j cannot support a central black hole when they have an 
isotropic dynamical structure, because stars at each binding en- 
ergy level will induce a r -1 ^ 2 cusp at small radii (Tremaine et 
al. 1994). 

5.3. The general case (0 < p < 1) 

For the general case, a direct calculation of the augmented 
density is not the best way to calculate the distribution func- 
tion. Instead, we use an approach based on the analysis of 
Ciotti (1996), who showed a convenient way to calculate the 
distribution function of a set of two-component models. If we 
invert the relation to r = r( v P s ,), and substitute it into the 
expression (|5}, we obtain 



¥ = o)Qi>,) = (1 + 



fia-t) 



(33) 



with t = f^*) as in expression d26l >. The relation J33l > links 
the stellar potential at a given position in the system to the 
total potential *P. We can also regard it more generally as the 
definition of a mapping function u, which maps the interval 
[0, T^] onto the positive real axis. It is a monotonically in- 
creasing function, and therefore the inverse function cj~ 1 exists. 
Transforming the Eddington formula J22I to an integration with 
respect to the stellar potential we obtain 



/(£) = 
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with £, 



1 (£) and Q a function defined by the integral 
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This expression is similar to the formula used by Ciotti (1996) 
to calculate the distribution function of his two-component 
Hernquist models. Whether this integration can be performed 



analytically depends on the complexity of both factors of the 
integrand. It appears that this integration can be performed an- 
alytically (only) for the models with y = 1, I and |. For these 
three models, both the augmented density p( v f,) and the map- 
ping function w( v F„) are rational functions. For y = 1 and 
y = |, the factor under the square root can be reduced to a cu- 
bic polynomial, while for y = | it can be reduced to a quadratic 
polynomial in ^T*. As a result, the function Q(x) and the distri- 
bution function can be expressed by means of elliptic integrals 
for the former two models, and completely in terms of elemen- 
tary functions for the latter model. More details can be found 
in Appendix B. 

For other values of y (including y = 2), an analytical eval- 
uation of the integral ( I34bi is not possible, and the distribution 
function has to be calculated numerically. For this goal, the ex- 
pression J34a> is not particularly useful, because it involves a 
differentiation of a numerically determined function. A more 
convenient formula for numerical integration can be obtained 
by using the alternative form fl23i of the Eddington equation. 
The second term in this expression vanishes for all y-models, 
because p(T') oc *P 4 at large radii. If we do the substitution 
W — > in this last equation, we obtain 



/(£) 



i r 6 * d r dp , &u 

~ VPJo dT^idWdT: 



dT, 
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Combined with a numerical solution of the equation £» = 
<xr'(£), this formula allows to evaluate the distribution func- 
tion numerically using standard quadrature techniques. 

In order to calculate the density of states function for y- 
models with a black hole, we can apply the same technique as 
for the calculation of the distribution function. If we rewrite 
equation ( 1241 as an integral with the stellar potential as the in- 
tegration variable, we obtain 



g(£)= 16 V2 



° 2 L 



,dr_ 



y/c4 y ¥ t )-oj(S t )d x V t . (36) 



This integral shows many similarities with the integral in the 
expression J34bl i. and the prospects to find an analytical solu- 
tion are very similar. One can demonstrate that the density of 
states function can be written in terms of elliptic integrals for 
y — 1 and y = |, and completely in terms of elementary func- 
tions for y = |. For all other values of y, including the Jaffe 
model, the density of states function cannot be expressed ana- 
lytically. 



5.4. Asymptotic expansions 

We can get a better insight into the behaviour of the distribution 
function and the DED of the y-models with a central black hole 
by studying the asymptotic behaviour. In the limit £ — > 0, i.e. 
in the outer regions of the system, a black hole has no effect 
on the dynamics of the y-models, as the asymptotic expansion 
of the potential *P — > r~' is independent of ju. The behaviour 
of the distribution function and the DED in the low binding- 
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Fig. 1. The distribution function (top row) and differential energy distribution (bottom row) for y-models without and with a 
central black hole. The value of y is indicated in each panel. The black curves correspond to the self-consistent models, the 
colored curves correspond to models with various black hole masses: = 0.001 (yellow), yi = 0.01 (green), fi = 0.1 (cyan) and 
H = 1 (red). 
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More interesting is the asymptotic behaviour of the distribu- 
tion function and the differential energy distribution in the high 
binding-energy limit. To find asymptotic expansions for the 
self-consistent y-models, we must separately consider the mod- 
els with a modest and those with a strong central density cusp. 
The self-consistent y-models with a modest density cusp y < 2 
have a finite potential well *To = (2 - yY l , and after some al- 
gebra one finds that in the limit £ — > ^o, 



V2 y(3-y) r (^) fcz 
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The self-consistent models with a central density cusp slope 
y > 2 have an infinitely deep potential well, and for these mod- 
els one finds in the limit £ — > oo 



/«(£) 



V2 J-y 6 -? 

y(3-y)(y-2)— £~, 



8tt 5 /2 



r(^) 



2y(3-y)rfe)r(|x) 



(41) 



(42) 



Finally, for the Jaffe model y = 2, the model that separates 
these two classes, the expansion in the limit £ — > 00 reads 



/.(£) 



1 



47r 5/2 



2 V6 o 
M(£) ~ e" s . 



(43) 
(44) 



When the y-models contain a black hole, the potential well is 
infinitely deep for all values of the cusp slope y. The asymp- 
totic behaviour of the distribution function and the differential 
energy distribution changes to 



/(£) 



N(S) 



3-7 r(y+l) 



2(2^/2 T(7-i) 
V^(3-y) rXy+l) 

8 r(7-|)' 



'£^ 3/2 , 



?- y s y ~ A . 



(45) 



(46) 



These expressions are valid for all models with a central black 
hole (0 < fi < 1) and for all values of y > \. For the limiting 
case y - \, the leading terms from equations ( 145 \ and ( 1461 
vanish and the asymptotic expansion for £ — > 00 reads 



V8 1 „, 
7V(S) ~4^ 5/2 £- 5 . 



(47) 
(48) 



5.5. Discussion 



In the top panels of figure^we plot the distribution function of 
various y-models with various black hole masses. A black hole 
drastically changes the behaviour of the distribution function, 
particularly for the models with a shallow cusp slope (y < 2) 
which have a finite stellar potential well. For these models, 
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the self-consistent distribution function is a strongly increas- 
ing function of energy, which diverges for S — > according 
to formula J39I . When a central black hole is present in these 
model, stars of all binding energies populate the galaxy. The 
behaviour of the distribution function in this new territory de- 
pends on the cusp slope, as prescribed by the asymptotic ex- 
pression J45i . For models with y < |, the distribution func- 
tion decreases as a function of binding energy in the neighbor- 
hood of the black hole. The distribution function of such mod- 
els thus has two regimes: it converges to zero both in the low 
and high binding energy limit and has a maximal value at bind- 
ing energies around the depth of the stellar potential well. The 
larger the black hole mass, the larger the value of the binding 
energy where the distribution function becomes maximal and 
the smoother the transition between the two regimes. For the 
y ' = I model, the distribution function becomes asymptotically 
fiat in the high binding energy limit. For models with a steeper 
cusp slope, the distribution function is a monotonically increas- 
ing function of binding energy, and the differences between a 
model without and with black hole becomes less pronounced. 
In particular for models with y > 2, the presence of a black hole 
does not drastically change the behaviour of the distribution 
function. Both without and with a central black hole, the distri- 
bution function is a monotonically increasing function of bind- 
ing energy, diverging in the high energy limit. Equations J41i 
and (145 1 show that the slope of the distribution function in the 
high energy limit changes from (6 - y)j (2y - 4) to y - | when 
a black hole is present. 

In the bottom row of figure^we plot the differential energy 
distribution N(S) for the same models as in the top row. For 
all models, M(S) converges to a finite value 8(3 - y)/5 in the 
low binding energies limit. Typically, N(S) is hardly affected 
for low binding energies, where it is a decreasing function of 
increasing black hole mass. Only at high binding energies, the 
effect of a black hole becomes visible, in particular for models 
with a shallow density cusp. The differential energy distribution 
N*(S) of the models without black hole suddenly drops to zero 
when £ approaches the depth of the potential well, in spite of 
the divergence of the distribution function [see equation ( I40H . 
In the presence of a black hole, where stars can populate or- 
bits with arbitrarily high binding energies, N(S) smoothly de- 
creases as £ r ~ 4 in the high energy limit. For y-models with a 
steeper cusp slope, in particular the models with an infinitely 
deep stellar potential well, the effect on the differential energy 
distribution is weaker. 

Although the effect of a black hole on the shape of the dis- 
tribution function is severe, it thus appears that the effect on the 
differential energy distribution is rather modest, even for the 
models with a finite stellar potential well. In the energy region 
that was off limits for the self-consistent models, N(S) assumes 
very low values. Although the addition of a black hole opens up 
a huge range for possible binding energies, the number of stars 
that actually populate these orbits is thus fairly small. This can 
be quantified by calculating the mean binding energy (£} of 
stars in the models. The mean binding energy is defined as 



-f 

Jo 



A more straightforward way to calculate the mean binding en- 
ergy is through the formula 

<£> = Iff ^ Iff ^ (r) " (50) 

Xoo 
p(r) [*F(r) - |cr 2 (r)] r 2 dr. 



(51) 



Using similar algebra as in section 4, we find 



4(5 - 2y) 



1 +jU 



4 - y 

3(2 -y) 



(52) 



Numerical integration of the equation ( I49l > for a selection of 
models gives identical results. For a Hernquist model (y = 1), 
the ratio (£)/(£)» equals 1 + fx, so the mean binding energy 
for a model with jj = 0.001 is only 0.1% larger than the mean 
binding energy of the self-consistent model. 

Another characteristic that can be used to quantify the im- 
portance of a central black hole on the energy distribution is the 
fraction 8 of stars on orbits with a binding energy larger than 
the potential well of the corresponding self-consistent model, 



N(8)d£. 



(53) 



For a y = 1 model with black hole masses fi = 0.001 and 
yu = 0.01, a numerical integration yields 9 = 1.63 x 10~ 4 and 
6 = 1.81 x 10~ 3 respectively. 

These numbers clearly demonstrate that the effect of a black 
hole on the global energy distribution is actually quite small, al- 
though the distribution function when represented as a function 
of the binding energy is affected in a very significant way. 

6. The LOSVD 

The LOSVD (also called line profile or velocity profile) 
<p p (R, V||) is the distribution of line-of-sight velocities at a given 
projected radius R. From an observer's point of view, this is 
definitely one of the most important kinematical quantities of 
a galaxy model, because a LOSVD contains all kinematic in- 
formation that can be obtained from a galaxy at a given line of 
sight and LOSVDs are in principle directly observable. 

6.1. Calculation of the LOSVDs 

The LOSVD can be found through the formula 



<P P ( R ^ v 'll) 



I(R) 



f 

Jr 



p(r)(/>(r,v\{)rdr 



(54) 



N(S)SdS. 



(49) 



where the function <p(r, v\\) represents the spatial LOSVD. In 
general, the spatial LOSVD <p{r, k, v\\) describes the distribu- 
tion of line-of-sight velocities at a position r in an arbitrary 
direction k. It can be found by integrating the distribution func- 
tion over the two velocity components perpendicular to k and 
normalizing the resulting distribution 

<P(r, k, v„) = ^y ) ff f( r > v ) (55) 




Fig. 2. The observed LOSVDs for y-models without and with a central black hole. The top panels are the LOSVDs at a projected 
radius R = 1, the middle panels correspond to R = 0.1 and the bottom panels are the LOSVDs at R — 0.01. The value of y is 
indicated each panel and the colour code is the same as in Figure^ 



For general anisotropic systems, these integrations are usually 
very cumbersome, and an analytical evaluation of the spatial 
LOSVD can only be obtained for a limited number of mod- 
els (e.g. Carollo et al. 1995; Baes & Dejonghe 2002b). For 
isotropic models, however, the spatial LOSVD is independent 
of the direction k and can be written as <p(r, v\\). The expres- 
sion d55i can then be transformed to 



<f>(r, V||) = 




where y - *T(r) - vjj/2. If we substitute the Eddington for- 
mula ( E2l into this equation, we obtain 

1 f v dp dY 

0(r,V||) = — (57) 

V27rp(r) Jo d*P ^Jy^ 

Comparing this expression with the Eddington formula ( 1221 . 
we see that the spatial LOSVD of an isotropic dynamical model 
can be calculated analytically if the same is true for the distribu- 
tion function. In fact, the spatial LOSVD is obtained almost au- 
tomatically during the calculation of the distribution function. 
This means that the spatial LOSVD of the self-consistent y- 
models can generally be expressed in terms of hypergeometric 
functions, and in terms of elementary functions for y — 2 ± l/n 
with n a natural number. The spatial LOSVD of the black hole 
dominated y-models can also be expressed in terms of hyper- 
geometric functions, and in terms of standard functions for 



integer and half-integer values of y. To calculate the spatial 
LOSVD for the y-models with a black hole, we use a simi- 
lar method as for the calculation of the distribution function. 
Changing the integration variable from *P to *P„ in equation d57l 
we find 

0(r, v N ) = — ± Q(y m ), (58) 

V27rp(r) 

with y» = aT l (y) and the function Q defined in equation J34bi . 
For the models with y = 1, y = | and y = |, both the inverse 
of the mapping function to and the function Q can be expressed 
analytically, such that the spatial LOSVD can be written com- 
pletely in analytical form (including elliptic integrals for the 
y = 1 and y = I models). For the other values of y, both of 
these operations must be done numerically. 

Summarizing, the calculation of the LOSVDs requires only 
one single quadrature for all self-consistent and black hole 
dominated models, and for the y-models with a black hole if 
y - 1 , y = |ory = | . For the remaining y-models containing 
a central black hole, the calculation of the LOSVDs requires a 
double quadrature. We calculated the LOSVDs with the Monte 
Carlo integration routines built into the SKIRT code (Baes et 
al. 2003b). 
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0.001 0.01 0.1 1 10 
projected radius r p 



0.0010.01 0.1 1 10 
projected radius r p 



0.0010.01 0.1 1 10 
projected radius r p 



0.001 0.01 0.1 1 10 
projected radius r p 



0.0010.01 0.1 1 10 
projected radius r p 



Fig. 3. The moments of the LOSVDs for y-models without and with a central black hole. The top row shows the projected 
dispersion profile (T P (R), the bottom row shows the Gauss-Hermite shape parameter h^iR). The value of y is indicated in each 
panel and the colour code is the same as in Figure^ 



6.2. Discussion 

In figure|21we plot a number of LOSVDs for a set of y-models 
for various values of the black hole mass fi. As expected, the in- 
fluence of a (realistic) black hole is negligible at large projected 
radii. At gradually smaller projected radii the influence of a 
black hole becomes more important. The signature of the pres- 
ence of the black hole is primarily a broadening of the LOS VD. 
The degree of the broadening depends on the black hole frac- 
tion and on the central slope. 

A quantitative investigation of the effects of a black hole on 
the LOSVDs is more easy when we study the moments of the 
LOSVDs. The first and second-order moments of the LOSVDs 
are nothing but the mean rotation and the projected dispersion. 
Rather than the true higher-order moments, one usually utilizes 
the coefficients h, from a Gauss-Hermite expansion to further 
characterize the shape of the LOSVDs (Gerhard 1993; van der 
Marel & Franx 1993). In figure[3]we plot the projected velocity 
dispersion profile and the I14 profile for a set of y-models with- 
out and with black holes of various masses (all odd moments 
are zero). 

The projected dispersion profile of models with a shallow 
central density slope are most sensitive to a black hole. As a re- 
sult of the presence of a black hole, the stars in the very centre 
of the galaxy will obtain large velocities, including velocities 
which are not possible in the finite central potential wells of the 
self-consistent models. This causes a dramatic increase in the 
velocity dispersion at small radii, resulting in the usual R~^ 2 
divergence. Note that this divergence is maintained in spite of 
the smoothing effect of the projection along the line of sight: 
the projected dispersion is a weighted integral of the intrinsic 
dispersion along the line of sight, such that the projected disper- 
sion at small projected radii contains a significant contribution 
from stars at large radii. 



The projected dispersion profile for models with a steep 
cusp slope is much less affected by a central black hole. Even 
without a black hole, these models already have steep stellar 
potential wells where stars can obtain arbitrarily high veloc- 
ities. The limit case is the (degenerate) model with y — > 3, 
where the self-consistent potential is already a point potential, 
and the addition of a black hole does not alter the kinematical 
structure of the galaxy at all. 

The effect of a black hole on the /14 shape parameter of 
the LOSVD is less straightforward to interpret. For the models 
with a shallow cusp, the central black hole does significantly 
affect the h\ profile, but there is no clear trend with increasing 
black hole mass. For the models with a steep density cusp, the 
effect of a black hole on the In profile is negligible. This trend, 
or rather the non-existence of a clear trend, contrasts with the 
signature of a dark matter halo on the observed kinematics at 
large radii. In elliptical galaxies, which lack a ubiquitous and 
straightforward tracer population such as the Hi gas in spiral 
galaxies, the stellar kinematics are one of the most important 
tracers for dark matter. The most obvious signature of dark mat- 
ter on the stellar kinematics is the behaviour of the projected 
velocity dispersion at large radii: dark matter causes the cr p 
profile to remain flat or to decrease more slowly than expected 
from the photometry alone (Saglia, Bertin & Stiavelli 1992; 
Saglia et al. 1993). This observational signature suffers from 
a mass-anisotropy degeneracy, however: a tangential stellar or- 
bital structure can cause a similar signature (Gerhard 1993). 
The key to discriminate between these possibilities is the be- 
haviour of the h\ profile. Carollo et al. (1995) studied in some 
detail the effect of a dark matter halo on the I14 shape parame- 
ter. They found that, both for isotropic and radially anisotropic 
models, the hn parameter increases at large projected radii due 
to the presence of a dark matter halo. On the other hand, a 
tangential orbital structure causes a negative I14 profile at large 
radii. The measurement of the projected dispersion and the h\ 
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profile at large radii would therefore in principle suffice to con- 
strain the dark matter content in ellipticals (Rix et al. 1997; 
Gerhard et al. 1998; Kronawitter et al. 2000), although system- 
atic effects caused by dust, non-sphericity etc. might further 
complicate this issue (Baes & Dejonghe 2001; Sanchis, Lokas 
& Mamon 2004). There is apparently no analogy for this situ- 
ation for the detection of black holes in galactic nuclei: as the 
signature of a central black hole on the profile does not fol- 
low a clear trend, the projected velocity dispersion is by far the 
most obvious observational signature by which a black hole can 
be detected in non-rotating stellar nuclei. 

7. Discussion and conclusion 

Galactic nuclei observed in real galaxies show a variety of 
structure, with central density cusps ranging from flat cores 
to steep cusps. The main goal of this paper was to construct 
simple dynamical models for spherical systems with a central 
black hole, reflecting the variety of central structure observed 
in real galactic nuclei. We have performed a detailed study of a 
two-parameter family of dynamical models based on the fam- 
ily of y-models introduced by Dehnen (1993) and Tremaine 
et al. (1994). The family contains as parameters the slope of 
the central density cusp y and the ratio // of the central black 
hole mass to the total mass of the system. By varying these 
parameters, we have been able to study a very wide range of 
models, going from weakly cusped models to very centrally 
concentrated models with an infinite stellar potential well, and 
ranging from self-consistent models without black hole to sys- 
tems where the dynamical structure is completely dominated 
by the central black hole. We have only considered models with 
an isotropic dynamical structure, which do not cover the range 
of orbital structure observed in real galactic nuclei. It is pos- 
sible to generalize the isotropic models presented in this pa- 
per to anisotropic models with a constant anisotropy or with 
an Osipkov-Merritt type distribution function (Osipkov 1979; 
Merritt 1985). For such models, the calculation of the most im- 
portant dynamical properties is not much more demanding than 
for isotropic models. The construction of general anisotropic 
dynamical models is much more complicated, however, and 
falls beyond the scope of this paper. 

For this two-parameter family of isotropic dynamical mod- 
els, we have calculated the most important kinematical quan- 
tities, such as the intrinsic and projected velocity dispersions, 
the total energy budget, the distribution function, the differen- 
tial energy distribution and the LOSVDs. Many of these quan- 
tities could be calculated completely analytically for all values 
of y and \x. For the models without black hole and the black 
hole dominated models, the distribution function, the differen- 
tial energy distribution and the spatial LOSVDs can be calcu- 
lated completely analytically. The same is true for three distinct 
y-models with a central black hole - these three models each 
form an analytical one-parameter family with the central black 
hole mass as an explicit parameter. Although these three mod- 
els differ only by the slope of the central density profile, their 
kinematical properties are very different. We therefore believe 
that the set of y-models discussed in this paper, and in particu- 
lar the three completely analytical models, answer the need for 



dynamical models which are on the one hand relatively sim- 
ple and on the other hand reflect the range of central structure 
observed in the nuclei of real galaxies. 

As most of the kinematics can be calculated completely an- 
alytically, the presented models allow to investigate the effect 
of a central black hole on the kinematics of galactic nuclei. 
In general, the effect of a black hole depends on the central 
density cusp (the parameter y). Models with a steep density 
cusp are least affected by the presence of a black hole. They 
already have a steep and infinitely deep potential well and a 
strong central density concentration, and the kinematical effect 
of a black hole is quite marginal. The distribution function re- 
mains a monotonically increasing function of binding energy; 
only the slope in the high-energy limit is altered by the presence 
of a black hole. Models with a shallow central density cusp on 
the other hand are more strongly affected by the presence of a 
black hole. They have a less concentrated density profile and 
therefore only a finite central potential well. The influence of 
a black hole on the distribution function can be important. For 
example, models with y < | have distribution functions which 
tend to zero both in the low and high binding energies regime. 
The differential energy distribution, which gives a better physi- 
cal insight in the orbital structure of the system, is less affected 
by the presence of a black hole. For realistic black hole masses, 
the shift in the mean binding energy of the stars or the fraction 
of stars on orbits with binding energies exceeding the depth 
of the stellar potential well is marginal. Nevertheless, the dif- 
ferent behaviour of the distribution function is important, as it 
bears direct consequences for the stability of the y-models with 
a black hole. A monotonically rising distribution function as a 
function of binding energy is a sufficient condition for stability 
against radial and non-radial perturbations. This means that the 
y-models with a black hole with y > | are stable. The mod- 
els with a shallower cusp slope have a decreasing distribution 
function in the high binding energy limit, and their stability is 
hence uncertain. This interesting issue can only be investigated 
in a detailed N-body or linear mode analysis studies. 

Finally, we can wonder whether these three models are 
unique or whether other simple dynamical models could eas- 
ily be found for which most of the kinematical properties can 
be expressed analytically. It is rather straightforward to con- 
structing models with a black hole in which the (intrinsic and/or 
projected) velocity dispersions can be expressed analytically, as 
the velocity dispersion is just a linear function of the black hole 
mass. The construction of models where more complicated dy- 
namical properties such as the distribution function and the spa- 
tial LOSVDs can be expressed analytically is more difficult. A 
direct integration of the Eddington equation does not seem the 
most obvious way to proceed [however, see Baes & Dejonghe 
(2004) for a case where this is doable]. A more promising path 
is the idea of mapping the total potential onto the stellar poten- 
tial, proposed by Ciotti (1996) and adapted by us in the present 
paper. Nevertheless, the conditions on the mapping function in 
order to allow an analytical evaluation of the distribution func- 
tion are quite stringent, and they probably do not apply for a 
large set of models. Out of the range of y-models, we found that 
only three different models satisfy these conditions. We also 
searched for other models in the large family of the so-called 
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{a,/5, y)-models (Zhao 1996), an extension of the y-models, but 
no other models satisfied the necessary conditions. 
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Appendix A: The function W y (r) 

The function W y (r) is defined in equation (|9jl as 



W y (r) 



3-7 r r 1_2y dr 



3-y r 

4* I 



(l+ r )7-2r' 



(A.1) 



where it is assumed that r is a real number and the parameter y 
assumes values between and 3. Formally, it can conveniently 
be expressed in terms of the incomplete Beta function or the 
hypergeometric function. It can also be written completely in 
terms of elementary functions. If y is not equal to 1, |, 
we obtain an algebraic expression 



,2or§, 



\2+j-2y 



7 4tt Lj\j)2 + j-2y [ \l + r) 



For the four special values of y, we can obtain a closed 
sion for W y (r) through direct integration of equation 
through application of the reduction formulae of hyper; 
ric functions. The results are 



1 (l+r 

Wi(r) = — In 

27T \ r 



3 (l + r 
W 3/2 (r) = - — In' 
Zn 



25 + 52r + 42r 2 + 12 r 3 
24tt(1 + rf ' 
3+22r + 30r 2 + 12r 3 



3 l+r\ 1 
W 2 (r) = —In — — 1 + 



1 l+r 
W5/2ir) = --ln — 



87rr(l +r) 3 
-4r- 18 r 2 - 12 r 3 
8nr 2 (l+r) 2 ' 
1 -2r + 6r 2 + 12 r 3 
247tr 3 (l+r) 



(A.2) 



expres- 
A~T1 or 

*eomet- 



(A.3) 
(A.4) 
(A.5) 
(A.6) 
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Appendix B: The function Q(x) for some selected 
models 

B.1. The y=l model 

For the Hernquist model we find a simple augmented density, 



POP*) = r- 



1 Wt 



In 1 - y t 
and for the mapping function we find 

fJ*¥* 



If we substitute the expressions 
tion J34bl . we obtain 



1 -¥„ 



(B.l) 



(B.2) 



B.ll and dB.21 into equa- 



Q(x) = 



— r 



2nyjl-n Jo V(^-*»)(l-^.) 3 



(B.3) 



where £ — 1 +///(1 -//)(l-x). The integrand in this expression is 
basically the combination of a rational function and the square 
root of a cubic polynomial. The function <3 (and hence the dis- 
tribution function) can therefore be expressed completely in 
terms of elliptic integrals. An explicit form requires a reduc- 
tion of this elliptic integral to a combination of the standard 
elliptic integrals. For the Hernquist model, one finds that the 
distribution can be written as 



/(£) = Ai + A 2 F(0, k) + A 3 E(0, k), 



(B.4) 



where A,- are algebraic functions of G and /j, F and E are 
Legendre's incomplete elliptic integrals of the first and second 
kind, and the arguments of the elliptic integrals are 



1 + 



= arcsin 
(l-M)(l-£) 2 
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For the calculation of the density of states, similar formulae 
apply. Practically, it is easiest to use a symbolic mathematical 
package such as Maple or Mathematica. These packages return 
an explicit form for most elliptic integrals, and the expressions 
can be converted immediately to Fortran or C code using inter- 
nal conversion routines. 



where the roots of the cubic polynomial under the square root 
are 
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Since the integrand of Q is the product of a rational function 
and the square root of a cubic polynomial, the integral can be 
written in terms of elliptic integrals. 



B.3. The 7 = § model 

The 7 = § model has an augmented density very similar to that 
of the 7 = | model, 
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The mapping function wCF») is very simple 
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and substitution of these expressions into equation ( I34bi gives 



Q(x) = 



l 
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(B.13) 

where ( - —(x + 4/yu). As the integrand in this expression 
is the product of a rational function and the square root of a 
quadratic polynomial, the function Q (and the distribution func- 
tion) can be expressed completely in terms of elementary func- 
tions. Explicit expressions for the distribution function, as well 
as for the density of states function, in term of elementary func- 
tion can be found in Baes & Dejonghe (2004). 



B.2. The 7 = 5 model 

For 7 = | we obtain 
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w( v F Ht ) = ¥ 
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Substitution of ( IB.7I and JB.8I into formula J34bl yields 
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